function [Betahat, Beta_se, Ehat, Sigmahat]=var_est(ym,xm)
% PURPOSE: VAR estimation 
% Estimates Beta for each of N equations through OLS
% Y_t = Y_{t-1} * B_1 + Y_{t-2} * B_2 + ... + Y_{t-nlags} * B_{nlags}
% 
% -----------------------------------------------------
% USAGE: [reg,u]=var_est(ym,xm);
% where: 
% ym = (T x N)   matrix of T observations of the N dependent variables
% xm = (T x nkk) matrix of T observations of the (nkk=N*nlags) independent variables
% NOTE: no dimension checks imposed
% -----------------------------------------------------
% RETURNS: 
% Betahat = parameters (column j contains the coefficients of eq(j) - constant&trend come first)
% Beta_se = standard errors of the parameter estimates (same ordering as Betahat) 
% Ehat    = (T x N) matrix of residuals.
% Sigmahat= estimate of the residual variance-covariance matrix
% tstat   = t-statistic of the parameters
% -----------------------------------------------------
% written by:
% Andrea Tamoni
% andrea.tamoni@phd.unibocconi.it
% � Dec 2007

N=size(ym,2);
[nt,nkk] = size(xm);

Betahat=inv(xm'*xm)*xm'*ym;          %nlags*N by N matrix of the parameters 
%betahatv=reshape(Betahat,nkk*nm,1); %nlags*N*N by 1 vector of the parameters

Ehat=ym-xm*Betahat;
Sigmahat=Ehat'*Ehat/(nt-nkk); %cov(Ehat)

VcvBetahat = kron(Sigmahat,inv(xm'*xm));
Beta_se = sqrt(diag(VcvBetahat));
Beta_se = reshape(Beta_se, nkk, N);

tstat = Betahat./Beta_se;